
// standard headers
#include <algorithm>
#include <cmath>
#include <iomanip>

// boost headers
#include <boost/format.hpp>

// component headers
#include <stlhelper.h>
#define DEBUG_LEVEL 0
#include <debug_func.h>
#include <tree_ops.h>
#include <matrix.h>

using boost::format;
using std::cerr;
using std::fill;
using std::vector;
using std::string;
using bcpp_matrix::t_dense_matrix_f;


namespace bcpp_tree_ops
{

//________________________________________________________________________________________

//	do_newhampshire_format_str

//________________________________________________________________________________________
void do_newhampshire_format_str::operator()(const bcpp_tree::t_node* root)
{
	debug_func_cerr(5);
	os << std::setprecision(5) << std::fixed;
	os << '(';
	for (unsigned i = 0; i < root->get_children_count()-1; ++i)
	{
		recurse(root->get_child(i));
		os << ",";
	}
	recurse(root->get_child(root->get_children_count()-1));
	os << ");";
}

void do_newhampshire_format_str::recurse(const bcpp_tree::t_node* node)
{
	assert(node);
	debug_func_cerr(5);
	unsigned cnt_children = node->get_children_count();
	// is leaf
//    format leaf_format("%1$s:%2$.5f");
//	format close_format("%1$.5f");

	if (!cnt_children)
	{
		os << node->get_name() << ":" << node->get_length();
//		os << node->get_name() << ":" << std::fixed << std::setprecision(5) << node->get_length();
		return;
	}
	os << '(';

	// handle all but last child
	for (unsigned i = 0; i < cnt_children - 1; ++i)
	{
		recurse(node->get_child(i));
		os << ',';
	}

	// handle last child
	recurse(node->get_child(cnt_children-1));

	os << "):" << node->get_length();
}







//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	do_collect_orthologs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//________________________________________________________________________________________

//	store_orthologs()


//		store orthologs in the various 1 to many sets

//________________________________________________________________________________________
namespace
{
std::string enode_type_names[] = {	"",							// eNONE
									"Orphan Tax1 Node",         // eTAX1
									"Orphan Tax2 Node",         // eTAX2
									"Orphan Tax1 Gene",         // eTAX1_GENE
									"Orphan Tax2 Gene",         // eTAX2_GENE
									"Ancestral Node",           // eTAX1_TAX2
									"Ortholog Tax1 Gene",       // eTAX1_GENE_ORTH
									"Ortholog Tax2 Gene",       // eTAX2_GENE_ORTH
									"Ortholog Tax1 Node",       // eTAX1_ORTH
									"Ortholog Tax2 Node",       // eTAX2_ORTH
									"Ortholog Node"};           // eTAX1_TAX2_ORTH
// for sorting vectors of genes by their first gene
struct t_gene_index
{
	std::string gene;
	unsigned	index;
	t_gene_index(const std::string& g, unsigned l) : gene(g), index(l){}
};

bool operator < (const t_gene_index& a, const t_gene_index& b)
{return a.gene < b.gene;}

};




struct cmp_gene_set
{
	bool operator () (const t_gene_set& a, const t_gene_set& b) const
	{	return a.genes[0] < b.genes[0];	}
};

void t_ortholog_data::sort_data()
{
	// each ortholog group ordered by name
	for (unsigned i = 0; i < gene_sets.size(); ++i)
		sort(gene_sets[i].genes.begin(), gene_sets[i].genes.end());

	// sort ortholog group by first name in each ortholog group
	sort_with_swaps(gene_sets.begin(), gene_sets.end(), cmp_gene_set());
}

void do_collect_orthologs::sort_data()
{
	// sort by gene names
	for (unsigned i =0; i < eORTH_TYPE_CNT; ++i)
		ortholog_sets[i].sort_data();
}

void do_collect_orthologs::store_genes(t_ortholog_data* store, float len)
{
	//cerr << "{";
	// save total count of genes
	store->cnt_tax1 +=tax1_genes.size();
	store->cnt_tax2 +=tax2_genes.size();

	// save ensg
	store->gene_sets.push_back(t_gene_set());
	for (unsigned i = 0; i < tax1_genes.size(); ++i)
		store->gene_sets.back().genes.push_back(tax1_genes[i]->get_name());
	for (unsigned i = 0; i < tax2_genes.size(); ++i)
		store->gene_sets.back().genes.push_back(tax2_genes[i]->get_name());
	store->gene_sets.back().length = len;

	// clear for next
	tax2_genes.clear();
	tax1_genes.clear();

	//cerr << "}";
}
void do_collect_orthologs::store_orphaned_paralogs()
{
	//cerr << "<";
#ifndef NDEBUG
	// make sure no orthologous genes
	for (unsigned i = 0; i < tax1_genes.size(); ++i)
		assert(get_node_type(tax1_genes[i])  == eTAX1_GENE);
	for (unsigned i = 0; i < tax2_genes.size(); ++i)
//		assert(node_types[tax2_genes[i]->get_UID()] == eTAX2_GENE);
		assert(get_node_type(tax2_genes[i])  == eTAX2_GENE);

#endif
	if (!tax1_genes.size() && !tax2_genes.size())
	{
		//cerr << ">";
		return;
	}

	t_ortholog_data* store = 0;


	// tax2 only
	if (!tax1_genes.size())
		store = &ortholog_sets[e_ORPH_TAX2];


	// tax1 only
	else
	{
		assert(tax2_genes.size() == 0);
		store = &ortholog_sets[e_ORPH_TAX1];
	}
	store_genes(store, -1);
	//cerr << ">";
}


void do_collect_orthologs::store_orthologs(bcpp_tree::t_tree::t_node_type* node)
{
	//cerr << "(";
#ifndef NDEBUG
	// make sure all orthologous genes
	for (unsigned i = 0; i < tax1_genes.size(); ++i)
	{

		if (get_node_type(tax1_genes[i]) != eTAX1_GENE_ORTH)
		{
			node->print(std::cerr);
			tax1_genes[i]->print(std::cerr);
		}
		assert(get_node_type(tax1_genes[i]) == eTAX1_GENE_ORTH);
	}
	for (unsigned i = 0; i < tax2_genes.size(); ++i)
	{
		if (get_node_type(tax2_genes[i]) != eTAX2_GENE_ORTH)
		{
			node->print(std::cerr);
			tax2_genes[i]->print(std::cerr);
		}
		assert(get_node_type(tax2_genes[i]) == eTAX2_GENE_ORTH);
	}
#endif
	assert(tax1_genes.size() && tax2_genes.size());

	t_ortholog_data* store = 0;

	if (tax1_genes.size() > 1)
	{
		if (tax2_genes.size() > 1)
			store = &ortholog_sets[e_ORTH_M_TO_M];
		else
			store = &ortholog_sets[e_ORTH_M_TO_1];
	}
	else
	{
		if (tax2_genes.size() > 1)
			store = &ortholog_sets[e_ORTH_1_TO_M];
		else
			store = &ortholog_sets[e_ORTH_1_TO_1];
	}

	store_genes(store, lengths_to_leaves[node->get_UID()]);
	//cerr << ")";
}


//________________________________________________________________________________________

//	operator()

//________________________________________________________________________________________
void do_collect_orthologs::recurse(bcpp_tree::t_tree::t_node_type* node)
{
	debug_func_cerr( 6 );
	eNODETYPE node_type = get_node_type(node);
	// take care of adjacent leaves first
	// ancestral node or ortholog node
	if (node_type == eTAX1_TAX2_ORTH || node_type == eTAX1_TAX2)
		store_orphaned_paralogs();

	// start from the bottom up
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i));

	// save list of orthologs in different 1:m, 1:1 etc. lists
	if (node_type == eTAX1_TAX2_ORTH)
		store_orthologs(node);
	else if (node_type == eTAX1_TAX2)
		store_orphaned_paralogs();

	// save actual genes
	else if (node_type == eTAX1_GENE_ORTH || node_type == eTAX1_GENE)
	{
		debug_cerr(5, "tax1 gene");
		tax1_genes.push_back(node);
	}
	else if (node_type == eTAX2_GENE_ORTH || node_type == eTAX2_GENE)
	{
		debug_cerr(5, "tax2 gene");
		tax2_genes.push_back(node);
	}

	else if (node_type == eTAX1)
	{
		debug_cerr(5, "tax1 orphan");
	}
	else if (node_type == eTAX2)
	{
		debug_cerr(5, "tax2 orphan");
	}
}




//________________________________________________________________________________________

//	do_set_node_types

//________________________________________________________________________________________
do_set_node_types::do_set_node_types(const bcpp_tree::t_tree& tree,
										const std::string& s1,
										const std::string& s2)
						 :
						tax1_matching_str(s1),
						tax2_matching_str(s2),
						regex_tax1(s1),
						regex_tax2(s2),
						node_types(tree.get_node_count(), eNONE)
{
	recurse(tree.get_root());
}
void do_set_node_types::recurse(bcpp_tree::t_tree::t_node_type* node)
{
	// start from bottom
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i));

	// leaves
	if (node->is_leaf())
	{
		if (regex_search(node->get_name(), regex_tax1, boost::match_default))
			set_node_type(*node, eTAX1_GENE);
		else if (regex_search(node->get_name(), regex_tax2, boost::match_default))
			set_node_type(*node, eTAX2_GENE);
		else
		{
			cerr	<<	"Error:\nUnknown gene:[" <<	node->get_name() << "]\n\n";
			cerr	<< node->get_name()
					<< " , gene 1 : " << tax1_matching_str
					<< " , gene 2 : " << tax2_matching_str
					<< "\n";
			set_node_type(*node, eNONE);
		}

		return;
	}


	// get list of leaf node types
	unsigned which_node_types_used[eTAX1_TAX2_ORTH + 1];
	fill(which_node_types_used,  which_node_types_used + eTAX1_TAX2_ORTH + 1,  0);
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		which_node_types_used[get_node_type(*node->get_child(i))]++;

	// divide into orthologs, tax1 and tax2 on the basis of leaf nodes

	// = Ancestral_Node
	//		if child includes Ancestral_Node or Ortholog_Node
	if (which_node_types_used[eTAX1_TAX2] || which_node_types_used[eTAX1_TAX2_ORTH])
		set_node_type(*node, eTAX1_TAX2);

	// = Orphan Tax1 Node
	//		if child includes only "Orphan Tax1 Gene" or "Orphan Tax1 Node"
	else if ( (which_node_types_used[eTAX1] + which_node_types_used[eTAX1_GENE] ) &&
			  !(which_node_types_used[eTAX2] + which_node_types_used[eTAX2_GENE] ))
		set_node_type(*node, eTAX1);

	// = Orphan Tax2 Node
	//		if child includes only "Orphan Tax2 Gene" or "Orphan Tax2 Node"
	else if (!(which_node_types_used[eTAX1] + which_node_types_used[eTAX1_GENE] ) &&
			 (which_node_types_used[eTAX2] + which_node_types_used[eTAX2_GENE] ))
		set_node_type(*node, eTAX2);

	// = Ortholog_Node
	//		if child includes both tax1 and tax2
	else
	{
		set_node_type(*node, eTAX1_TAX2_ORTH);
		tag_ortholog_decendents(node);
	}
}
void do_set_node_types::tag_ortholog_decendents(bcpp_tree::t_tree::t_node_type* node)
{
	// tag orthologous duplication nodes as lineage specific duplications
	// rather than loses

	if (node_types[node->get_UID()] == eTAX1)
		node_types[node->get_UID()]=eTAX1_ORTH;
	else if (node_types[node->get_UID()] == eTAX2)
		node_types[node->get_UID()]=eTAX2_ORTH;
	else if (node_types[node->get_UID()] == eTAX1_GENE)
		node_types[node->get_UID()]=eTAX1_GENE_ORTH;
	else if (node_types[node->get_UID()] == eTAX2_GENE)
		node_types[node->get_UID()]=eTAX2_GENE_ORTH;
	assert (node_types[node->get_UID()] != eNONE);

	// recurse
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		tag_ortholog_decendents(node->get_child(i));
}






//________________________________________________________________________________________

//	do_calculate_branch_lengths

//________________________________________________________________________________________
do_calculate_branch_lengths::do_calculate_branch_lengths(const bcpp_tree::t_tree& tree)
	:lengths_to_leaves(tree.get_node_count(), 0.0),
	max_lengths_to_leaves(tree.get_node_count(), 0.0),
	 lengths_to_root(tree.get_node_count(), 0.0),
	 max_len_to_leaf(0.0)
{
	// set lengths from root
	recurse_len_to_root(tree.get_root());

	// set lengths from leaf
	float min_len_to_leaves = 0.0;
	float max_len_to_leaves = 0.0;
	recurse_len_to_leaves(tree.get_root(), min_len_to_leaves, max_len_to_leaves);
}

void do_calculate_branch_lengths::recurse_len_to_root(bcpp_tree::t_tree::t_node_type* node)
{
	float len_to_root = node->get_length();
	if (node->get_parent())
		len_to_root += lengths_to_root[node->get_parent()->get_UID()];

	lengths_to_root[node->get_UID()] = len_to_root;

	// recurse each child
	if (node->is_leaf())
		// save leaf length
		max_len_to_leaf = std::max(max_len_to_leaf, len_to_root);
	else
		for (unsigned i = 0; i < node->get_children_count(); ++i)
			recurse_len_to_root(node->get_child(i));
}

void do_calculate_branch_lengths::recurse_len_to_leaves(bcpp_tree::t_tree::t_node_type* node,
						float& parent_len_to_leaf_min, float& parent_len_to_leaf_max)
{
//	lengths_to_leaves[node->get_UID()] =
//									max_len_to_leaf - lengths_to_root[node->get_UID()];
//	// recurse each child
//	for (unsigned i = 0; i < node->get_children_count(); ++i)
//		recurse_len_to_leaves(node->get_child(i));

	if (node->is_leaf())
	{
		// length to leaves for itself is zero
		lengths_to_leaves[node->get_UID()] =
		max_lengths_to_leaves[node->get_UID()] = 0.0;

		// return length to leaf for its parent
		parent_len_to_leaf_min =
		parent_len_to_leaf_max =	node->get_length();
		return;
	}

	assert(node->get_children_count());
	parent_len_to_leaf_min = 10000.0;
	parent_len_to_leaf_max = 0.0;
	for (unsigned i = 0; i < node->get_children_count(); ++i)
	{
		// get min and max length to leaf via each child
		float leaf_via_child_min;
		float leaf_via_child_max;
		recurse_len_to_leaves(node->get_child(i), leaf_via_child_min, leaf_via_child_max);
		if (leaf_via_child_min < parent_len_to_leaf_min)
			parent_len_to_leaf_min = leaf_via_child_min;
		if (leaf_via_child_max > parent_len_to_leaf_max)
			parent_len_to_leaf_max = leaf_via_child_max;
	}
	// save length to leaf
	lengths_to_leaves[node->get_UID()] = parent_len_to_leaf_min;
	max_lengths_to_leaves[node->get_UID()] = parent_len_to_leaf_max;

	// return length to leaf for its parent
	parent_len_to_leaf_min += node->get_length();
	parent_len_to_leaf_max += node->get_length();
}



//________________________________________________________________________________________

//	op_print_tree

//________________________________________________________________________________________

op_print_tree::op_print_tree(	const bcpp_tree::t_tree&	tree,
								std::ostream&	 	 			os_,
								const std::vector<float>&		lengths_to_leaves_,
								const std::vector<float>&		lengths_to_root_,
								const std::vector<eNODETYPE>&	node_types_,
								bool							nodes_only):
				lengths_to_leaves(lengths_to_leaves_), lengths_to_root(lengths_to_root_),
				node_types(node_types_),
				os(os_), curr_level(0)
{
	if (nodes_only)
		recurse_nodes_only(tree.get_root());
	else
		// print tree
		recurse(tree.get_root());
}

void op_print_tree::recurse(bcpp_tree::t_tree::t_node_type* node)
{
	eNODETYPE node_type = node_types.size() ? node_types[node->get_UID()] : eNONE;
	os << format("%1$s--[%2$.5f / %3$.5f / %4$.5f]-->%5$s  {%6$s}\n")
		% string(curr_level * 5, ' ')
		% node->get_length()
		% lengths_to_root[node->get_UID()]
		% lengths_to_leaves[node->get_UID()]
		% node->get_name()
		% enode_type_names[node_type];
	++curr_level;
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i));
	--curr_level;
}



void op_print_tree::recurse_nodes_only(bcpp_tree::t_tree::t_node_type* node)
{
	eNODETYPE node_type = node_types.size() ? node_types[node->get_UID()] : eNONE;
	os << format("%1$.5f\t%2$.5f\t%3$s\t%4$s\n")
			% lengths_to_root[node->get_UID()]
			% lengths_to_leaves[node->get_UID()]
			% node->get_name()
			% enode_type_names[node_type];
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse_nodes_only(node->get_child(i));
}




//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	op_print_tree_elaborate

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
void op_print_tree_elaborate::print() const
{
	debug_func_cout( 5 );
	os.precision(precision);
	os << std::fixed;
	tree.for_each(*this);
}

void op_print_tree_elaborate::print(const bcpp_tree::t_tree::t_node_type& curr_node) const
{
	debug_func_cout( 5 );
	os.precision(precision);
	os << std::fixed;
	recurse(curr_node);
}
void op_print_tree_elaborate::recurse(const bcpp_tree::t_tree::t_node_type& curr_node) const
{
	if (curr_node.is_leaf())
	{
		print_horizontal_lines();
#ifndef NDEBUG
		if (node_types.size() <= curr_node.get_UID() && node_types.size())
		{
			cerr << "curr_node.get_UID() = " << curr_node.get_UID() << "\n";
		}
		assert(node_types.size() > curr_node.get_UID() || node_types.size() == 0);
#endif
		eNODETYPE node_type = node_types.size() ? node_types[curr_node.get_UID()] : eNONE;
		os << curr_node.get_name();
		if (node_type != eNONE)
			os << " {" << enode_type_names[node_type] << "}\n";
		else
			os << "\n";

		return;
	}
	else
	{
		print_horizontal_lines();
		levels.push_back(1);

		// first child
		if (curr_node.get_children_count())
		{
			assert(curr_node.get_child(0));
			recurse(*curr_node.get_child(0));
		}

		// middle children
		for (unsigned i = 1; i < curr_node.get_children_count(); ++i)
		{
			print_vertical_lines();
			// last child, print node distance
			//if (i == curr_node.get_children_count() - 1)
				os << lengths_to_leaves[curr_node.get_UID()];
			os << std::endl;
			unsigned last_level = levels.back();
			// last child, print with vertical bar
			if (i == curr_node.get_children_count() - 1)
				levels.back() = 0;
			print_vertical_lines();
			recurse(*curr_node.get_child(i));
			levels.back() = last_level;
		}

		levels.pop_back();
	}
}




















//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	do_sum_sqrs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


void do_sum_sqrs::operator()(bcpp_tree::t_tree::t_node_type* root)
{
	debug_func_cout(5);
	std::deque<unsigned> parent_branch;
	recurse(root, parent_branch);
	avg_percent_stddev	= cnt_non_null_data == 1 ? 0.0 :
						100.0 * sqrt(sum_sqrs / (cnt_non_null_data * 2.0 - 2.0));
//	APSD =sqrt ( SSQ / (N-2) )    x 100
}



void do_sum_sqrs::recurse(bcpp_tree::t_tree::t_node_type* node, std::deque<unsigned>& parent_branch)
{
	debug_func_cout(5);
	assert(node);

	// save index in parent's list if leaf
	if (node->is_leaf())
	{
		parent_branch.push_back(node->get_leaf_UID());
		return;
	}

	assert( node->get_children_count() == 2);

	// create lists of indices for all descendents of each child
	std::deque<std::deque<unsigned> >	branches(node->get_children_count());
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i), branches[i]);

	// length from node to leaf
	float node_dist = lengths_to_leaves[node->get_UID()] * 2;

	assert(branches.size());
	for (unsigned i = 0; i + 1 < branches.size(); ++i)
	{
		for (unsigned j = i + 1; j < branches.size(); ++j)
		{
			std::deque<unsigned>& indices1 = branches[i];
			std::deque<unsigned>& indices2 = branches[j];
			for (unsigned ii = 0; ii < indices1.size(); ++ii)
			{
				for (unsigned jj = 0; jj < indices2.size(); ++jj)
				{
					assert(indices1[ii] < matrix.size());
					assert(indices2[jj] < matrix.size());
					float D = matrix[indices1[ii]][indices2[jj]];
					// ignore null
					if (D <= 0 || D > 5000.0	)
						continue;
					sum_sqrs += 2 * (D - node_dist) * (D - node_dist) / (D *D);
					++cnt_non_null_data;
				}
			}
		}
	}

	// now move all into the branch of my parent
	for (unsigned i = 0; i < branches.size(); ++i)
		parent_branch.insert(parent_branch.begin(), branches[i].begin(), branches[i].end());
	return;
}
void CHECK_LINKS(bcpp_tree::t_tree::t_node_type& curr_node, std::ostream& os)
{
	static unsigned indent = 0;
	os << string(indent, ' ') << curr_node.get_name() << std::endl;
	indent += 2;
	if (!curr_node.is_leaf())
	{
		assert(curr_node.get_child(0));
		assert(curr_node.get_child(1));
		CHECK_LINKS(*curr_node.get_child(0), os);
		CHECK_LINKS(*curr_node.get_child(1), os);
	}
	indent -= 2;
}






//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	do_sum_sqrs2

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
void do_sum_sqrs2::operator()(bcpp_tree::t_tree::t_node_type* root)
{
	debug_func_cout(5);
	branch_lengths.assign(root->get_tree()->get_leaf_count(), 0.0);
	branch_differences.assign(root->get_tree()->get_leaf_count(), 0.0);
	std::deque<unsigned> parent_branch;
	recurse(root, parent_branch);
	avg_percent_stddev	= cnt_non_null_data == 1 ? 0.0 :
						100.0 * sqrt(sum_sqrs / (cnt_non_null_data * 2.0 - 2.0));
//	APSD =sqrt ( SSQ / (N-2) )    x 100
}



void do_sum_sqrs2::recurse(bcpp_tree::t_tree::t_node_type* node, std::deque<unsigned>& parent_branch)
{
//	debug_func_cout(1);
	assert(node);

	// save index in parent's list if leaf
	if (node->is_leaf())
	{
		assert((int)branch_lengths.size() > node->get_leaf_UID());
		parent_branch.push_back(node->get_leaf_UID());
		branch_lengths[node->get_leaf_UID()] += node->get_length();
		return;
	}

	assert( node->get_children_count() == 2);

	// create lists of indices for all descendents of each child
	std::deque<std::deque<unsigned> >	branches(node->get_children_count());

	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i), branches[i]);
	float node_dist = lengths_to_leaves[node->get_UID()] * 2;

	assert(branches.size());
	for (unsigned i = 0; i + 1 < branches.size(); ++i)
	{
		for (unsigned j = i + 1; j < branches.size(); ++j)
		{
			std::deque<unsigned>& indices1 = branches[i];
			std::deque<unsigned>& indices2 = branches[j];
			for (unsigned ii = 0; ii < indices1.size(); ++ii)
			{
				for (unsigned jj = 0; jj < indices2.size(); ++jj)
				{
					unsigned gene_1 = indices1[ii];
					unsigned gene_2 = indices2[jj];
					assert(fabs(node_dist - (branch_lengths[gene_2] +
											 branch_lengths[gene_1])) < 0.001);;

					assert(gene_1 < matrix.size());
					assert(gene_2 < matrix.size());
					float D = matrix[gene_1][gene_2];
					// ignore null
					if (D <= 0.0 || D > 5000.0	)
						continue;
					if (node_dist > 5000.0	)
						continue;
					float diff = 2 * (D - node_dist) * (D - node_dist) / (D *D);
					sum_sqrs += diff;
					branch_differences[gene_1] += diff;
					branch_differences[gene_2] += diff;
					++cnt_non_null_data;

				}
			}
		}
	}

	// now move all into the branch of my parent
	for (unsigned i = 0; i < branches.size(); ++i)
		parent_branch.insert(parent_branch.begin(), branches[i].begin(), branches[i].end());
	for (unsigned i = 0; i < parent_branch.size(); ++i)
		branch_lengths[parent_branch[i]] += node->get_length();
	return;
}









//________________________________________________________________________________________

//	op_extract_branch_leaves

//________________________________________________________________________________________
void op_extract_branch_leaves::recurse(const bcpp_tree::t_tree::t_node_type* node)
{
	if (node->is_leaf())
	{
		leaves.push_back(node);
		return;
	}
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i));
}


//________________________________________________________________________________________

//	t_species_tree_map

//________________________________________________________________________________________
void t_species_tree_map::recurse(bcpp_tree::t_node* node, std::deque<unsigned>& child_UIDs, unsigned& taxa_count)
{
	unsigned UID = node->get_UID();
	child_UIDs.push_back(UID);

	//
	// If leaf, then in-species duplications counted?
	//

	if (node->is_leaf())
	{
		// missing taxa = maximum???
		cnt_missing_taxa[UID][UID] = 0;
		// duplications = 1???
		cnt_duplications[UID][UID] = 1;
		// same-species paralogs
		parent_species_UID[UID][UID] = UID;
		cnt_child_taxa_at_node[UID] = 1;
		++taxa_count;
		return;
	}


	//
	// Recurse to get children data
	//

	// all children nodes (both taxa and intermediate nodes)
	std::deque<unsigned> left_UIDs, rite_UIDs;
	// count of all taxa in children
	unsigned child_taxa_count = 0;
	assert(node->get_children_count() == 2);
	recurse(node->get_child(0),  left_UIDs, child_taxa_count);
	recurse(node->get_child(1),  rite_UIDs, child_taxa_count);


	// increment leaf count by the number of child taxa
	taxa_count += cnt_child_taxa_at_node[UID] = child_taxa_count;


	//
	// children vs children relationships
	//

	for (deq_iter ii = left_UIDs.begin(); ii != left_UIDs.end(); ++ii)
	{
		for (deq_iter jj = rite_UIDs.begin(); jj != rite_UIDs.end(); ++jj)
		{
			// missing genes if taxa skipped
			cnt_missing_taxa[*ii][*jj] =
			cnt_missing_taxa[*jj][*ii] = child_taxa_count
											- cnt_child_taxa_at_node[*ii]
											- cnt_child_taxa_at_node[*jj];
			// not orphan
			cnt_duplications[*ii][*jj] =
			cnt_duplications[*jj][*ii] = 0;
			// when children join, resulting node is me
			parent_species_UID[*ii][*jj] =
			parent_species_UID[*jj][*ii] = UID;
		}
	}


	// return all child node UIDs in child_UIDs
	child_UIDs.insert(child_UIDs.begin(),  left_UIDs.begin(),  left_UIDs.end());
	child_UIDs.insert(child_UIDs.begin(),  rite_UIDs.begin(),  rite_UIDs.end());

	//
	// My relationships vs children and myself
	// Parent_species_UID is always me
	// Obviously duplications (cnt_duplications defaults to 1)
	// 	and hence species missing, but how many?
	//

	for (deq_iter ii = child_UIDs.begin(); ii != child_UIDs.end(); ++ii)
	{
		cnt_missing_taxa[*ii][UID] =
		cnt_missing_taxa[UID][*ii] = child_taxa_count - cnt_child_taxa_at_node[*ii];
		parent_species_UID[*ii][UID] =
		parent_species_UID[UID][*ii] = UID;

	}

}


//________________________________________________________________________________________

//	op_set_labels_to_children_names

//________________________________________________________________________________________
void op_set_labels_to_children_names::recurse(bcpp_tree::t_node* node, std::string& label)
{
	if (node->is_leaf())
	{
		node->set_label(node->get_name());
		if (label.length())
			label += '-' + node->get_label();
		else
			label = node->get_label();
		return;
	}
	string child_label;
	recurse(node->get_child(0), child_label);
	recurse(node->get_child(1), child_label);

	node->set_label(child_label);
	if (label.length())
		label += '-' + child_label;
	else
		label = child_label;
}


//________________________________________________________________________________________

//	op_tree_names_to_short_names

//________________________________________________________________________________________
void op_tree_names_to_short_names::recurse(bcpp_tree::t_tree::t_node_type& node)
{
	if (node.is_leaf())
	{
		string long_name = node.get_name();
		string short_name = (boost::format("XYZX%06d") % current_index++).str();
		node.set_name(short_name);
		long_name_to_short_name[long_name] = short_name;
	}
	else
		for (unsigned i = 0; i < node.get_children_count(); ++i)
			recurse(*(node.get_child(i)));
}







//________________________________________________________________________________________

//	op_dump_tree_data

//________________________________________________________________________________________

op_dump_tree_data::op_dump_tree_data(	const bcpp_tree::t_tree&	tree_,
										std::ostream&	 	 		os_):
				os(os_), curr_level(0), tree(tree_)
{
}

void op_dump_tree_data::dump()
{
	tree.assert_consistency();

	// print tree
	os << "tree_name = " << tree.get_name()				<< "\n";
	os << "number of leaves = " << tree.get_leaf_count()		<< "\n";
	os << "number of nodes  = " << tree.get_node_count()		<< "\n";
	os << "root pointer= " << reinterpret_cast<int*>(tree.get_root())	<< "\n";
	recurse(tree.get_root());
}

void op_dump_tree_data::recurse(bcpp_tree::t_tree::t_node_type* node)
{
	os << format(
		"%|1$|--"
		"%|2$| "				// leaf/node
		"\"%|3$|"					// name
		"%|4$|\", "				// label
		"[%|5$g|], "			// length
		"UID=%|6$d|"			// UID
		"(leafUID=%|7$d|), "	// leaf UID
		"children=%|8$d| "		// children
		"addr=%|9$d|<-%|10$d| ")	// pointer
		% string(curr_level * 5, ' ')
		% (node->is_leaf() ? "leaf " : "node ")
		% node->get_name()
		% (node->get_label().length() ? "(" + node->get_label() + ")" : "")
		% node->get_length()
		% node->get_UID()
		% node->get_leaf_UID()
		% node->get_children_count()
		% long(node->get_parent())
		% long(node);
	if (node->get_children_count())
	{
		os << "("	<< (long) (node->get_child(0));
		for (unsigned i = 1; i < node->get_children_count(); ++i)
			os <<  "," << (long) (node->get_child(i));
		os << ")";
	}
	os << "\n";
	
	++curr_level;
	for (unsigned i = 0; i < node->get_children_count(); ++i)
		recurse(node->get_child(i));
	--curr_level;
}

};

